% StiffeningConcentrationProfile_3D_PEG.m
% Eden Ford
% 02 August 2025

% Code to model discretized diffusion-reaction concentration profile into a
% cylindrical hydrogel


% Please excuse the lack of efficiency in this code. I am not a computer
% scientist.

% 1 second time steps
% update C_inf every hour

close all
clear all
clc

C_XT_max = 12.32; % mM, for plotting C_XT

%% Variable Definitions
% C_A = fixed azide concentration (fixed to gel)
% C_B = fixed BCN concentration (fixed to gel)
% C_P0 = free PEG concentration (diffusing)
% C_PF = fixed PEG concentration (no longer diffusing)
% C_P = total PEG concentration (free + fixed; both influence free PEG 
%   diffusion)
% C_F = functional group on PEG for reaction (assume that if all BCN and
%   azides react, and all PEGs have the same number of reactions per
%   PEG, only this concentration of functional groups is 'available'
%   for reaction
% C_X = new SPACC linkage concentration (fixed to gel)
% C_XT = sum of C_X over the total incubations (shows total increase in
%   linkage density
% u = C_F/C_P; constant incorporated into governing equations to represent
%   the concentration of functional groups 'available' for reaction

%% General Constants and Variables
% All in cm, s, mmol

D_water = 4.81e-07; % cm^2/s; diffusivity in water for PEG
D = 2.80e-07; % cm^2/s; effective diffusivity for PEG
k = 0.57; % 1/(M*s) [=] cm3/(mmol*s); SPAAC rate constant at 37C
V = 0.25; % mL [=] cm3; solution volume

R = 4.6/10; % cm; hydrogel radius
H = 2.0/10; % cm; hydrogel height

delta_r = 0.005; % cm; 50 um step size
delta_z = 0.005; % cm; 50 um step size
delta_t = 1; % seconds; 1 second intervals for time steps

r_list = [0 : delta_r : R]; % cm; range of r locations
z_list = [0 : delta_z : H]; % cm; range of z locations
[r,z] = meshgrid(r_list,z_list);
N = length(r_list);
M = length(z_list);

%% Stiffening Step 1A (PEG-4BCN): Initial Conditions

writerPEG = VideoWriter('PEG_Video');
writerPEG.FrameRate = 2;
open(writerPEG)

writerXL = VideoWriter('XL_Video');
writerXL.FrameRate = 2;
open(writerXL)

disp('Stiffening Step 1A: PEG-4BCN');

C_A0 = 0.77/1000; % mmol/cm3 [=] mM/1000; fixed azide functional groups throughout gel
C_i_BCN = 4.68/1000; % mmol/cm3 [=] mM/1000; initial BCN concentration in solution (250 mL)
C_i = C_i_BCN/4; % mmol/cm3 [=] mM/1000; initial PEG-4BCN concentration in solution (250 mL)

A_max = C_A0; % for plotting

% NOTE: treat PEG as though only one BCN will react and the other 3 will be
% fixed to gel
u = 1; % 1 'available' BCN for reaction per PEG
t_f = 0;
disp(['Timepoint: ' num2str(t_f) ' hours'])

EmptyMatrix = zeros(size(r));
C_A = EmptyMatrix+C_A0; % mmol/cm3; fixed azide concentration
C_P0 = EmptyMatrix; % mmol/cm3; free PEG-4BCN concentration
C_PF = EmptyMatrix; % mmol/cm3; newly fixed PEG-4BCN concentration
C_X = EmptyMatrix; % mmol/cm3; new XL concentration
C_XT = EmptyMatrix; % mmol/cm3; total new XL concentration over course of entire stiffening

C_P = C_P0+C_PF; % mmol/cm3; total new PEG concentration

C_P0_plot = 1000*C_P0; % mM; free PEG-4BCN concentration
C_PF_plot = 1000*C_PF; % mM; fixed PEG-4BCN concentration
C_P_plot = 1000*C_P;
C_X_plot = 1000*C_X; % mM; new XL concentration (this stiffening step)
C_XT_plot = 1000*(C_XT+C_X); % mM; total new XL concentration over course of entire stiffening

% Heatmap Plot
fig = figure();
fig.Position(3:4) = [1300 650];
hold on
    % free PEG concentration (C_P0)
    ax1 = subplot(3,3,3);
    Plot1 = surf(r,z,C_P0_plot);
    set(Plot1,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(autumn);
    colormap(ax1,cmap)
    title(['Solution Concentration: ' num2str(C_i*1000) ' mM PEG'],'FontSize',14)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'free PEG concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 1.755])
    txt = ['Free PEG'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
    
    % fixed PEG Concentration (C_PF)
    ax2 = subplot(3,3,6);
    Plot2 = surf(r,z,C_PF_plot);
    set(Plot2,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(spring);
    colormap(ax2,cmap)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'fixed PEG concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 1.155])
    txt = ['Fixed PEG'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
    
    % total PEG (C_P)
    ax3 = subplot(3,3,9);
    Plot3 = surf(r,z,C_P_plot);
    set(Plot3,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(copper);
    colormap(ax3,cmap)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'total PEG concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 2.91])
    txt = ['Total PEG'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
    
    % stiffened linkage concentration (C_X)
    ax4 = subplot(3,3,2);
    Plot4 = surf(r,z,C_X_plot);
    set(Plot4,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(gray);
    colormap(ax4,cmap)
    title('Initial Conditions','FontSize',14)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'new linkage concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([-0.05 2.31])
    txt = ['New Linkages (Within Incubation)'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
hold off
frame = getframe(gcf);
writeVideo(writerPEG,frame)

% Heatmap Plot
fig = figure();
fig.Position(3:4) = [800 500];
hold on
    % total linkage concentration (C_XT)
    Plot = surf(r,z,C_XT_plot);
    set(Plot,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    title('Initial Conditions','FontSize',16)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'total stiffened linkage concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 C_XT_max])
    txt = ['Total New Covalent Linkages Over the Course of Stiffening'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
hold off
frame = getframe(gcf);
writeVideo(writerXL,frame)

%% Stiffening Step 1A (PEG-4BCN)
C_inf = C_i; % mmol/cm3 [=] mM/1000;  initial PEG-4BCN concentration in solution (250 mL)
disp(['Solution Concentration: ' num2str(C_inf*1000) ' mM PEG'])

PlotTime = [1 : 1 : 12]; % in hours
% Cycle through PlotTime
for b = 1:length(PlotTime)
    time = PlotTime(b);
    t0 = t_f + delta_t; % initial timepoint (in seconds)
    t_f = time*60*60; % plot timepoint in seconds
    t_list = [t0 : delta_t : t_f]; % t steps
    disp(['Timepoint: ' num2str(time) ' hours'])
    
    % Calculate concentrations at each time step
    for a = 1:length(t_list)
        C_A_NEW = C_A-k*delta_t*C_A.*C_P0*u; % new fixed azide concentration
        C_X_NEW = C_X+k*delta_t*C_A.*C_P0*u; % new XL concentration
        C_PF_NEW = C_PF+k*delta_t*C_A.*C_P0; % new fixed PEG concentration
        
        % Calculate new free PEG concentration
        C_P0_NEW = EmptyMatrix;
        
        % BC at r = R (i = N)
        C_P0_NEW(:,end) = C_inf;
        % BC at z = H (j = M)
        C_P0_NEW(end,:) = C_inf;
        % BC at r = 0 (i = 1), z = 0 (j = 1)
        C_P0_NEW(1,1) = C_P0(1,1) - ((7*D*delta_t*C_P(1,1))/(delta_r^2)+(7*D*delta_t*C_P(1,1))/(2*delta_z^2)+k*delta_t*C_A(1,1)*C_P0(1,1)) + ((8*D*delta_t)/(delta_r^2))*C_P(1,2) - ((D*delta_t)/(delta_r^2))*C_P(1,3) + ((4*D*delta_t)/(delta_z^2))*C_P(2,1) - ((D*delta_t)/(2*delta_z^2))*C_P(3,1);
        % 'BC' (axis of symmetry) at r = 0 (i = 1)
        for j = 2:M-1
            C_P0_NEW(j,1) = C_P0(j,1) - ((7*D*delta_t*C_P(j,1))/(delta_r^2)+(2*D*delta_t*C_P(j,1))/(delta_z^2)+k*delta_t*C_A(j,1)*C_P0(j,1)) + ((8*D*delta_t)/(delta_r^2))*C_P(j,2) - ((D*delta_t)/(delta_r^2))*C_P(j,3) + ((D*delta_t)/(delta_z^2))*C_P(j+1,1) + ((D*delta_t)/(delta_z^2))*C_P(j-1,1);
        end
        % BC at z = 0 (j = 1)
        for i = 2:N-1
            r_i = (i-1)*delta_r;
            C_P0_NEW(1,i) = C_P0(1,i) - ((2*D*delta_t*C_P(1,i))/(delta_r^2)+(7*D*delta_t*C_P(1,i))/(2*delta_z^2)+k*delta_t*C_A(1,i)*C_P0(1,i)) + ((4*D*delta_t)/(delta_z^2))*C_P(2,i) - ((D*delta_t)/(2*delta_z^2))*C_P(3,i) + ((D*delta_t)/(delta_r^2) + (D*delta_t)/(2*r_i*delta_r))*C_P(1,i+1) + ((D*delta_t)/(delta_r^2) - (D*delta_t)/(2*r_i*delta_r))*C_P(1,i-1);
        end
        % internal nodes
        for j = 2:M-1
            for i = 2:N-1
                C_P0_NEW(j,i) = C_P0(j,i) - ((2*D*delta_t*C_P(j,i))/(delta_r^2)+(2*D*delta_t*C_P(j,i))/(delta_z^2)+k*delta_t*C_A(j,i)*C_P0(j,i)) + ((D*delta_t)/(delta_r^2) + (D*delta_t)/(2*r_i*delta_r))*C_P(j,i+1) + ((D*delta_t)/(delta_r^2) - (D*delta_t)/(2*r_i*delta_r))*C_P(j,i-1) + ((D*delta_t)/(delta_z^2))*C_P(j+1,i) + ((D*delta_t)/(delta_z^2))*C_P(j-1,i);
            end
        end
        
        % Update new concentrations
        C_A = C_A_NEW; % mmol/cm3; fixed azide concentration
        C_P0 = C_P0_NEW; % mmol/cm3; free PEG-4BCN concentration
        C_PF = C_PF_NEW; % mmol/cm3; newly fixed PEG-4BCN concentration
        C_X = C_X_NEW; % mmol/cm3; new XL concentration
        C_P = C_P0+C_PF; % mmol/cm3; total new PEG concentration
    end
        
    % Update bulk solution concentration
    C_gel_slice = zeros(M-1,N-1);
    for i = 1:N-1
        for j = 1:M-1
            C_gel_slice(j,i) = delta_r*delta_z*((C_P(j,i)+C_P(j,i+1)+C_P(j+1,i)+C_P(j+1,i+1))/4); % mmol/cm3
            % this is the PEG (fixed + free) concentration along a
            % slice of the gel: at each box of delta r and delta z
            % concentration of free PEG + reacted PEG
        end
    end
    C_gel_slice_sum = sum(sum(C_gel_slice)); % equivalent to C_P*R*H (mmol/cm)
    C_gel_full = C_gel_slice_sum*pi*R; % remainder of volume (mmol)
    C_inf = (C_i*V - C_gel_full)/V; % mmol/cm3; new C_inf value for next hour
    
    disp(['New Solution Concentration: ' num2str(C_inf*1000) ' mM PEG'])
    
    % Prepare to plot
    C_P0_plot = 1000*C_P0; % mM; free PEG-4BCN concentration
    C_PF_plot = 1000*C_PF; % mM; fixed PEG-4BCN concentration
    C_P_plot = 1000*C_P; % mM; total PEG-4BCN concentration
    C_X_plot = 1000*C_X; % mM; new XL concentration (this stiffening step)
    C_XT_plot = 1000*(C_XT+C_X); % mM; total new XL concentration over course of entire stiffening
    
    % Heatmap Plot
    fig = figure();
    fig.Position(3:4) = [1300 650];
    hold on
        % free PEG concentration (C_P0)
        ax1 = subplot(3,3,3);
        Plot1 = surf(r,z,C_P0_plot);
        set(Plot1,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(autumn);
        colormap(ax1,cmap)
        title(['Solution Concentration: ' num2str(C_inf*1000) ' mM PEG'],'FontSize',14)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'free PEG concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 1.755])
        txt = ['Free PEG'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off

        % fixed PEG Concentration (C_PF)
        ax2 = subplot(3,3,6);
        Plot2 = surf(r,z,C_PF_plot);
        set(Plot2,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(spring);
        colormap(ax2,cmap)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'fixed PEG concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 1.155])
        txt = ['Fixed PEG'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off

        % total PEG (C_P)
        ax3 = subplot(3,3,9);
        Plot3 = surf(r,z,C_P_plot);
        set(Plot3,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(copper);
        colormap(ax3,cmap)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'total PEG concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 2.91])
        txt = ['Total PEG'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off

        % stiffened linkage concentration (C_X)
        ax4 = subplot(3,3,2);
        Plot4 = surf(r,z,C_X_plot);
        set(Plot4,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(gray);
        colormap(ax4,cmap)
        title([num2str(time) '-Hour Timepoint'],'FontSize',14)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'new linkage concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([-0.05 2.31])
        txt = ['New Linkages (Within Incubation)'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off
    hold off
    frame = getframe(gcf);
    writeVideo(writerPEG,frame)

    % Heatmap Plot
    fig = figure();
    fig.Position(3:4) = [800 500];
    hold on
        % total linkage concentration (C_XT)
        Plot = surf(r,z,C_XT_plot);
        set(Plot,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        title([num2str(time) '-Hour Timepoint'],'FontSize',14)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'total stiffened linkage concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 C_XT_max])
        txt = ['Total New Covalent Linkages Over the Course of Stiffening'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off
    hold off
    frame = getframe(gcf);
    writeVideo(writerXL,frame)
end
C_B = 3*C_X; % mmol/cm3; fixed BCN concentration
C_XT = C_X+C_XT; % update total linkages over course of stiffening

P0_min = min(min(C_P0))*1000
P0_max = max(max(C_P0))*1000
PF_min = min(min(C_PF))*1000
PF_max = max(max(C_PF))*1000
P_min = min(min(C_P))*1000
P_max = max(max(C_P))*1000
X_min = min(min(C_X))*1000
X_max = max(max(C_X))*1000
XT_min = min(min(C_XT))*1000
XT_max = max(max(C_XT))*1000
%pause(120)

%% Stiffening Step 1B (PEG-4Azide): Initial Conditions

disp(' ');
disp('Stiffening Step 1B: PEG-4Azide');

C_B0 = C_B; % mmol/cm3 [=] mM/1000; fixed BCN functional groups throughout gel
C_A0 = C_A; % mmol/cm3 [=] mM/1000; fixed azide functional groups throughout gel
% Disregard any free PEG remaining in the gel

C_i_Az = 7.02/1000; % mmol/cm3 [=] mM/1000; initial azide concentration in solution (250 mL)
C_i = C_i_Az/4; % mmol/cm3 [=] mM/1000; initial PEG-4Az concentration in solution (250 mL)

B_max = max(max(C_B0)); % for plotting

% NOTE: treat PEG as though only 2 azides will react and the other 2 will
% be fixed to gel
u = 2; % 2 'available' azide for reaction per PEG
t_f = 0;
disp(['Timepoint: ' num2str(t_f) ' hours'])

EmptyMatrix = zeros(size(r));
C_B = EmptyMatrix+C_B0; % mmol/cm3; fixed BCN concentration
C_P0 = EmptyMatrix; % mmol/cm3; free PEG-4Az concentration
C_PF = EmptyMatrix; % mmol/cm3; newly fixed PEG-4Az concentration
C_X = EmptyMatrix; % mmol/cm3; new XL concentration

C_P = C_P0+C_PF; % mmol/cm3; total new PEG concentration

C_P0_plot = 1000*C_P0; % mM; free PEG-4Az concentration
C_PF_plot = 1000*C_PF; % mM; fixed PEG-4Az concentration
C_P_plot = 1000*C_P;
C_X_plot = 1000*C_X; % mM; new XL concentration (this stiffening step)
C_XT_plot = 1000*(C_XT+C_X); % mM; total new XL concentration over course of entire stiffening

% Heatmap Plot
fig = figure();
fig.Position(3:4) = [1300 650];
hold on
    % free PEG concentration (C_P0)
    ax1 = subplot(3,3,3);
    Plot1 = surf(r,z,C_P0_plot);
    set(Plot1,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(autumn);
    colormap(ax1,cmap)
    title(['Solution Concentration: ' num2str(C_i*1000) ' mM PEG'],'FontSize',14)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'free PEG concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 1.755])
    txt = ['Free PEG'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
    
    % fixed PEG Concentration (C_PF)
    ax2 = subplot(3,3,6);
    Plot2 = surf(r,z,C_PF_plot);
    set(Plot2,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(spring);
    colormap(ax2,cmap)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'fixed PEG concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 1.155])
    txt = ['Fixed PEG'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
    
    % total PEG (C_P)
    ax3 = subplot(3,3,9);
    Plot3 = surf(r,z,C_P_plot);
    set(Plot3,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(copper);
    colormap(ax3,cmap)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'total PEG concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 2.91])
    txt = ['Total PEG'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
    
    % stiffened linkage concentration (C_X)
    ax4 = subplot(3,3,2);
    Plot4 = surf(r,z,C_X_plot);
    set(Plot4,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(gray);
    colormap(ax4,cmap)
    title('Initial Conditions','FontSize',14)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'new linkage concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([-0.05 2.31])
    txt = ['New Linkages (Within Incubation)'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
hold off
frame = getframe(gcf);
writeVideo(writerPEG,frame)

% Heatmap Plot
fig = figure();
fig.Position(3:4) = [800 500];
hold on
    % total linkage concentration (C_XT)
    Plot = surf(r,z,C_XT_plot);
    set(Plot,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    title('Initial Conditions','FontSize',16)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'total stiffened linkage concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 C_XT_max])
    txt = ['Total New Covalent Linkages Over the Course of Stiffening'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
hold off
frame = getframe(gcf);
writeVideo(writerXL,frame)

%% Stiffening Step 1B (PEG-4Azide)
C_inf = C_i; % mmol/cm3 [=] mM/1000;  initial PEG-4Az concentration in solution (250 mL)
disp(['Solution Concentration: ' num2str(C_inf*1000) ' mM PEG'])

PlotTime = [1 : 1 : 12]; % in hours
% Cycle through PlotTime
for b = 1:length(PlotTime)
    time = PlotTime(b);
    t0 = t_f + delta_t; % initial timepoint (in seconds)
    t_f = time*60*60; % plot timepoint in seconds
    t_list = [t0 : delta_t : t_f]; % t steps
    disp(['Timepoint: ' num2str(time) ' hours'])
    
    % Calculate concentrations at each time step
    for a = 1:length(t_list)
        C_B_NEW = C_B-k*delta_t*C_B.*C_P0*u; % new fixed BCN concentration
        C_X_NEW = C_X+k*delta_t*C_B.*C_P0*u; % new XL concentration
        C_PF_NEW = C_PF+k*delta_t*C_B.*C_P0; % new fixed PEG concentration
        
        % Calculate new free PEG concentration
        C_P0_NEW = EmptyMatrix;
        
        % BC at r = R (i = N)
        C_P0_NEW(:,end) = C_inf;
        % BC at z = H (j = M)
        C_P0_NEW(end,:) = C_inf;
        % BC at r = 0 (i = 1), z = 0 (j = 1)
        C_P0_NEW(1,1) = C_P0(1,1) - ((7*D*delta_t*C_P(1,1))/(delta_r^2)+(7*D*delta_t*C_P(1,1))/(2*delta_z^2)+k*delta_t*C_B(1,1)*C_P0(1,1)) + ((8*D*delta_t)/(delta_r^2))*C_P(1,2) - ((D*delta_t)/(delta_r^2))*C_P(1,3) + ((4*D*delta_t)/(delta_z^2))*C_P(2,1) - ((D*delta_t)/(2*delta_z^2))*C_P(3,1);
        % 'BC' (axis of symmetry) at r = 0 (i = 1)
        for j = 2:M-1
            C_P0_NEW(j,1) = C_P0(j,1) - ((7*D*delta_t*C_P(j,1))/(delta_r^2)+(2*D*delta_t*C_P(j,1))/(delta_z^2)+k*delta_t*C_B(j,1)*C_P0(j,1)) + ((8*D*delta_t)/(delta_r^2))*C_P(j,2) - ((D*delta_t)/(delta_r^2))*C_P(j,3) + ((D*delta_t)/(delta_z^2))*C_P(j+1,1) + ((D*delta_t)/(delta_z^2))*C_P(j-1,1);
        end
        % BC at z = 0 (j = 1)
        for i = 2:N-1
            r_i = (i-1)*delta_r;
            C_P0_NEW(1,i) = C_P0(1,i) - ((2*D*delta_t*C_P(1,i))/(delta_r^2)+(7*D*delta_t*C_P(1,i))/(2*delta_z^2)+k*delta_t*C_B(1,i)*C_P0(1,i)) + ((4*D*delta_t)/(delta_z^2))*C_P(2,i) - ((D*delta_t)/(2*delta_z^2))*C_P(3,i) + ((D*delta_t)/(delta_r^2) + (D*delta_t)/(2*r_i*delta_r))*C_P(1,i+1) + ((D*delta_t)/(delta_r^2) - (D*delta_t)/(2*r_i*delta_r))*C_P(1,i-1);
        end
        % internal nodes
        for j = 2:M-1
            for i = 2:N-1
                C_P0_NEW(j,i) = C_P0(j,i) - ((2*D*delta_t*C_P(j,i))/(delta_r^2)+(2*D*delta_t*C_P(j,i))/(delta_z^2)+k*delta_t*C_B(j,i)*C_P0(j,i)) + ((D*delta_t)/(delta_r^2) + (D*delta_t)/(2*r_i*delta_r))*C_P(j,i+1) + ((D*delta_t)/(delta_r^2) - (D*delta_t)/(2*r_i*delta_r))*C_P(j,i-1) + ((D*delta_t)/(delta_z^2))*C_P(j+1,i) + ((D*delta_t)/(delta_z^2))*C_P(j-1,i);
            end
        end
        
        % Update new concentrations
        C_B = C_B_NEW; % mmol/cm3; fixed BCN concentration
        C_P0 = C_P0_NEW; % mmol/cm3; free PEG-4Az concentration
        C_PF = C_PF_NEW; % mmol/cm3; newly fixed PEG-4Az concentration
        C_X = C_X_NEW; % mmol/cm3; new XL concentration
        C_P = C_P0+C_PF; % mmol/cm3; total new PEG concentration
    end
        
    % Update bulk solution concentration
    C_gel_slice = zeros(M-1,N-1);
    for i = 1:N-1
        for j = 1:M-1
            C_gel_slice(j,i) = delta_r*delta_z*((C_P(j,i)+C_P(j,i+1)+C_P(j+1,i)+C_P(j+1,i+1))/4); % mmol/cm3
            % this is the PEG (fixed + free) concentration along a
            % slice of the gel: at each box of delta r and delta z
            % concentration of free PEG + reacted PEG
        end
    end
    C_gel_slice_sum = sum(sum(C_gel_slice)); % equivalent to C_P*R*H (mmol/cm)
    C_gel_full = C_gel_slice_sum*pi*R; % remainder of volume (mmol)
    C_inf = (C_i*V - C_gel_full)/V; % mmol/cm3; new C_inf value for next hour

    disp(['New Solution Concentration: ' num2str(C_inf*1000) ' mM PEG'])
    
    % Prepare to plot
    C_P0_plot = 1000*C_P0; % mM; free PEG-4BCN concentration
    C_PF_plot = 1000*C_PF; % mM; fixed PEG-4BCN concentration
    C_P_plot = 1000*C_P; % mM; total PEG-4BCN concentration
    C_X_plot = 1000*C_X; % mM; new XL concentration (this stiffening step)
    C_XT_plot = 1000*(C_XT+C_X); % mM; total new XL concentration over course of entire stiffening
    
    % Heatmap Plot
    fig = figure();
    fig.Position(3:4) = [1300 650];
    hold on
        % free PEG concentration (C_P0)
        ax1 = subplot(3,3,3);
        Plot1 = surf(r,z,C_P0_plot);
        set(Plot1,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(autumn);
        colormap(ax1,cmap)
        title(['Solution Concentration: ' num2str(C_inf*1000) ' mM PEG'],'FontSize',14)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'free PEG concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 1.755])
        txt = ['Free PEG'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off

        % fixed PEG Concentration (C_PF)
        ax2 = subplot(3,3,6);
        Plot2 = surf(r,z,C_PF_plot);
        set(Plot2,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(spring);
        colormap(ax2,cmap)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'fixed PEG concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 1.155])
        txt = ['Fixed PEG'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off

        % total PEG (C_P)
        ax3 = subplot(3,3,9);
        Plot3 = surf(r,z,C_P_plot);
        set(Plot3,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(copper);
        colormap(ax3,cmap)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'total PEG concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 2.91])
        txt = ['Total PEG'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off

        % stiffened linkage concentration (C_X)
        ax4 = subplot(3,3,2);
        Plot4 = surf(r,z,C_X_plot);
        set(Plot4,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(gray);
        colormap(ax4,cmap)
        title([num2str(time) '-Hour Timepoint'],'FontSize',14)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'new linkage concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([-0.05 2.31])
        txt = ['New Linkages (Within Incubation)'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off
    hold off
    frame = getframe(gcf);
    writeVideo(writerPEG,frame)

    % Heatmap Plot
    fig = figure();
    fig.Position(3:4) = [800 500];
    hold on
        % total linkage concentration (C_XT)
        Plot = surf(r,z,C_XT_plot);
        set(Plot,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        title([num2str(time) '-Hour Timepoint'],'FontSize',14)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'total stiffened linkage concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 C_XT_max])
        txt = ['Total New Covalent Linkages Over the Course of Stiffening'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off
    hold off
    frame = getframe(gcf);
    writeVideo(writerXL,frame)
end
C_A = C_X+C_A0; % mmol/cm3; fixed azide concentration (new + leftover from last cycle)
C_XT = C_X+C_XT; % update total linkages over course of stiffening

P0_min = min(min(C_P0))*1000
P0_max = max(max(C_P0))*1000
PF_min = min(min(C_PF))*1000
PF_max = max(max(C_PF))*1000
P_min = min(min(C_P))*1000
P_max = max(max(C_P))*1000
X_min = min(min(C_X))*1000
X_max = max(max(C_X))*1000
XT_min = min(min(C_XT))*1000
XT_max = max(max(C_XT))*1000
%pause(120)

%% Stiffening Step 2A (PEG-4BCN): Initial Conditions

disp(' ');
disp('Stiffening Step 2A: PEG-4BCN');

C_B0 = C_B; % mmol/cm3 [=] mM/1000; fixed BCN functional groups throughout gel
C_A0 = C_A; % mmol/cm3 [=] mM/1000; fixed azide functional groups throughout gel
% Disregard any free PEG remaining in the gel

C_i_BCN = 7.02/1000; % mmol/cm3 [=] mM/1000; initial BCN concentration in solution (250 mL)
C_i = C_i_BCN/4; % mmol/cm3 [=] mM/1000; initial PEG-4BCN concentration in solution (250 mL)

A_max = max(max(C_A0)); % for plotting

% NOTE: treat PEG as though only 2 BCNs will react and the other 2 will be
% fixed to gel
u = 2; % 2 'available' BCN for reaction per PEG
t_f = 0;
disp(['Timepoint: ' num2str(t_f) ' hours'])

EmptyMatrix = zeros(size(r));
C_A = EmptyMatrix+C_A0; % mmol/cm3; fixed azide concentration
C_P0 = EmptyMatrix; % mmol/cm3; free PEG-4BCN concentration
C_PF = EmptyMatrix; % mmol/cm3; newly fixed PEG-4BCN concentration
C_X = EmptyMatrix; % mmol/cm3; new XL concentration

C_P = C_P0+C_PF; % mmol/cm3; total new PEG concentration

C_P0_plot = 1000*C_P0; % mM; free PEG-4Az concentration
C_PF_plot = 1000*C_PF; % mM; fixed PEG-4Az concentration
C_P_plot = 1000*C_P;
C_X_plot = 1000*C_X; % mM; new XL concentration (this stiffening step)
C_XT_plot = 1000*(C_XT+C_X); % mM; total new XL concentration over course of entire stiffening

% Heatmap Plot
fig = figure();
fig.Position(3:4) = [1300 650];
hold on
    % free PEG concentration (C_P0)
    ax1 = subplot(3,3,3);
    Plot1 = surf(r,z,C_P0_plot);
    set(Plot1,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(autumn);
    colormap(ax1,cmap)
    title(['Solution Concentration: ' num2str(C_i*1000) ' mM PEG'],'FontSize',14)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'free PEG concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 1.755])
    txt = ['Free PEG'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
    
    % fixed PEG Concentration (C_PF)
    ax2 = subplot(3,3,6);
    Plot2 = surf(r,z,C_PF_plot);
    set(Plot2,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(spring);
    colormap(ax2,cmap)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'fixed PEG concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 1.155])
    txt = ['Fixed PEG'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
    
    % total PEG (C_P)
    ax3 = subplot(3,3,9);
    Plot3 = surf(r,z,C_P_plot);
    set(Plot3,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(copper);
    colormap(ax3,cmap)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'total PEG concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 2.91])
    txt = ['Total PEG'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
    
    % stiffened linkage concentration (C_X)
    ax4 = subplot(3,3,2);
    Plot4 = surf(r,z,C_X_plot);
    set(Plot4,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(gray);
    colormap(ax4,cmap)
    title('Initial Conditions','FontSize',14)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'new linkage concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([-0.05 2.31])
    txt = ['New Linkages (Within Incubation)'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
hold off
frame = getframe(gcf);
writeVideo(writerPEG,frame)

% Heatmap Plot
fig = figure();
fig.Position(3:4) = [800 500];
hold on
    % total linkage concentration (C_XT)
    Plot = surf(r,z,C_XT_plot);
    set(Plot,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    title('Initial Conditions','FontSize',16)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'total stiffened linkage concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 C_XT_max])
    txt = ['Total New Covalent Linkages Over the Course of Stiffening'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
hold off
frame = getframe(gcf);
writeVideo(writerXL,frame)

%% Stiffening Step 2A (PEG-4BCN)
C_inf = C_i; % mmol/cm3 [=] mM/1000;  initial PEG-4BCN concentration in solution (250 mL)
disp(['Solution Concentration: ' num2str(C_inf*1000) ' mM PEG'])

PlotTime = [1 : 1 : 12]; % in hours
% Cycle through PlotTime
for b = 1:length(PlotTime)
    time = PlotTime(b);
    t0 = t_f + delta_t; % initial timepoint (in seconds)
    t_f = time*60*60; % plot timepoint in seconds
    t_list = [t0 : delta_t : t_f]; % t steps
    disp(['Timepoint: ' num2str(time) ' hours'])
    
    % Calculate concentrations at each time step
    for a = 1:length(t_list)
        C_A_NEW = C_A-k*delta_t*C_A.*C_P0*u; % new fixed azide concentration
        C_X_NEW = C_X+k*delta_t*C_A.*C_P0*u; % new XL concentration
        C_PF_NEW = C_PF+k*delta_t*C_A.*C_P0; % new fixed PEG concentration
        
        % Calculate new free PEG concentration
        C_P0_NEW = EmptyMatrix;
        
        % BC at r = R (i = N)
        C_P0_NEW(:,end) = C_inf;
        % BC at z = H (j = M)
        C_P0_NEW(end,:) = C_inf;
        % BC at r = 0 (i = 1), z = 0 (j = 1)
        C_P0_NEW(1,1) = C_P0(1,1) - ((7*D*delta_t*C_P(1,1))/(delta_r^2)+(7*D*delta_t*C_P(1,1))/(2*delta_z^2)+k*delta_t*C_A(1,1)*C_P0(1,1)) + ((8*D*delta_t)/(delta_r^2))*C_P(1,2) - ((D*delta_t)/(delta_r^2))*C_P(1,3) + ((4*D*delta_t)/(delta_z^2))*C_P(2,1) - ((D*delta_t)/(2*delta_z^2))*C_P(3,1);
        % 'BC' (axis of symmetry) at r = 0 (i = 1)
        for j = 2:M-1
            C_P0_NEW(j,1) = C_P0(j,1) - ((7*D*delta_t*C_P(j,1))/(delta_r^2)+(2*D*delta_t*C_P(j,1))/(delta_z^2)+k*delta_t*C_A(j,1)*C_P0(j,1)) + ((8*D*delta_t)/(delta_r^2))*C_P(j,2) - ((D*delta_t)/(delta_r^2))*C_P(j,3) + ((D*delta_t)/(delta_z^2))*C_P(j+1,1) + ((D*delta_t)/(delta_z^2))*C_P(j-1,1);
        end
        % BC at z = 0 (j = 1)
        for i = 2:N-1
            r_i = (i-1)*delta_r;
            C_P0_NEW(1,i) = C_P0(1,i) - ((2*D*delta_t*C_P(1,i))/(delta_r^2)+(7*D*delta_t*C_P(1,i))/(2*delta_z^2)+k*delta_t*C_A(1,i)*C_P0(1,i)) + ((4*D*delta_t)/(delta_z^2))*C_P(2,i) - ((D*delta_t)/(2*delta_z^2))*C_P(3,i) + ((D*delta_t)/(delta_r^2) + (D*delta_t)/(2*r_i*delta_r))*C_P(1,i+1) + ((D*delta_t)/(delta_r^2) - (D*delta_t)/(2*r_i*delta_r))*C_P(1,i-1);
        end
        % internal nodes
        for j = 2:M-1
            for i = 2:N-1
                C_P0_NEW(j,i) = C_P0(j,i) - ((2*D*delta_t*C_P(j,i))/(delta_r^2)+(2*D*delta_t*C_P(j,i))/(delta_z^2)+k*delta_t*C_A(j,i)*C_P0(j,i)) + ((D*delta_t)/(delta_r^2) + (D*delta_t)/(2*r_i*delta_r))*C_P(j,i+1) + ((D*delta_t)/(delta_r^2) - (D*delta_t)/(2*r_i*delta_r))*C_P(j,i-1) + ((D*delta_t)/(delta_z^2))*C_P(j+1,i) + ((D*delta_t)/(delta_z^2))*C_P(j-1,i);
            end
        end
        
        % Update new concentrations
        C_A = C_A_NEW; % mmol/cm3; fixed azide concentration
        C_P0 = C_P0_NEW; % mmol/cm3; free PEG-4BCN concentration
        C_PF = C_PF_NEW; % mmol/cm3; newly fixed PEG-4BCN concentration
        C_X = C_X_NEW; % mmol/cm3; new XL concentration
        C_P = C_P0+C_PF; % mmol/cm3; total new PEG concentration
    end
        
    % Update bulk solution concentration
    C_gel_slice = zeros(M-1,N-1);
    for i = 1:N-1
        for j = 1:M-1
            C_gel_slice(j,i) = delta_r*delta_z*((C_P(j,i)+C_P(j,i+1)+C_P(j+1,i)+C_P(j+1,i+1))/4); % mmol/cm3
            % this is the PEG (fixed + free) concentration along a
            % slice of the gel: at each box of delta r and delta z
            % concentration of free PEG + reacted PEG
        end
    end
    C_gel_slice_sum = sum(sum(C_gel_slice)); % equivalent to C_P*R*H (mmol/cm)
    C_gel_full = C_gel_slice_sum*pi*R; % remainder of volume (mmol)
    C_inf = (C_i*V - C_gel_full)/V; % mmol/cm3; new C_inf value for next hour
    
    disp(['New Solution Concentration: ' num2str(C_inf*1000) ' mM PEG'])
    
    % Prepare to plot
    C_P0_plot = 1000*C_P0; % mM; free PEG-4BCN concentration
    C_PF_plot = 1000*C_PF; % mM; fixed PEG-4BCN concentration
    C_P_plot = 1000*C_P; % mM; total PEG-4BCN concentration
    C_X_plot = 1000*C_X; % mM; new XL concentration (this stiffening step)
    C_XT_plot = 1000*(C_XT+C_X); % mM; total new XL concentration over course of entire stiffening
    
    % Heatmap Plot
    fig = figure();
    fig.Position(3:4) = [1300 650];
    hold on
        % free PEG concentration (C_P0)
        ax1 = subplot(3,3,3);
        Plot1 = surf(r,z,C_P0_plot);
        set(Plot1,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(autumn);
        colormap(ax1,cmap)
        title(['Solution Concentration: ' num2str(C_inf*1000) ' mM PEG'],'FontSize',14)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'free PEG concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 1.755])
        txt = ['Free PEG'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off

        % fixed PEG Concentration (C_PF)
        ax2 = subplot(3,3,6);
        Plot2 = surf(r,z,C_PF_plot);
        set(Plot2,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(spring);
        colormap(ax2,cmap)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'fixed PEG concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 1.155])
        txt = ['Fixed PEG'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off

        % total PEG (C_P)
        ax3 = subplot(3,3,9);
        Plot3 = surf(r,z,C_P_plot);
        set(Plot3,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(copper);
        colormap(ax3,cmap)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'total PEG concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 2.91])
        txt = ['Total PEG'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off

        % stiffened linkage concentration (C_X)
        ax4 = subplot(3,3,2);
        Plot4 = surf(r,z,C_X_plot);
        set(Plot4,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(gray);
        colormap(ax4,cmap)
        title([num2str(time) '-Hour Timepoint'],'FontSize',14)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'new linkage concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([-0.05 2.31])
        txt = ['New Linkages (Within Incubation)'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off
    hold off
    frame = getframe(gcf);
    writeVideo(writerPEG,frame)

    % Heatmap Plot
    fig = figure();
    fig.Position(3:4) = [800 500];
    hold on
        % total linkage concentration (C_XT)
        Plot = surf(r,z,C_XT_plot);
        set(Plot,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        title([num2str(time) '-Hour Timepoint'],'FontSize',14)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'total stiffened linkage concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 C_XT_max])
        txt = ['Total New Covalent Linkages Over the Course of Stiffening'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off
    hold off
    frame = getframe(gcf);
    writeVideo(writerXL,frame)
end
C_B = C_X+C_B0; % mmol/cm3; fixed BCN concentration
C_XT = C_X+C_XT; % update total linkages over course of stiffening

P0_min = min(min(C_P0))*1000
P0_max = max(max(C_P0))*1000
PF_min = min(min(C_PF))*1000
PF_max = max(max(C_PF))*1000
P_min = min(min(C_P))*1000
P_max = max(max(C_P))*1000
X_min = min(min(C_X))*1000
X_max = max(max(C_X))*1000
XT_min = min(min(C_XT))*1000
XT_max = max(max(C_XT))*1000
%pause(120)

%% Stiffening Step 2B (PEG-4Azide): Initial Conditions

disp(' ');
disp('Stiffening Step 2B: PEG-4Azide');

C_B0 = C_B; % mmol/cm3 [=] mM/1000; fixed BCN functional groups throughout gel
C_A0 = C_A; % mmol/cm3 [=] mM/1000; fixed azide functional groups throughout gel
% Disregard any free PEG remaining in the gel

C_i_Az = 7.02/1000; % mmol/cm3 [=] mM/1000; initial azide concentration in solution (250 mL)
C_i = C_i_Az/4; % mmol/cm3 [=] mM/1000; initial PEG-4Az concentration in solution (250 mL)

B_max = max(max(C_B0)); % for plotting

% NOTE: treat PEG as though only 2 azides will react and the other 2 will
% be fixed to gel
u = 2; % 2 'available' azide for reaction per PEG
t_f = 0;
disp(['Timepoint: ' num2str(t_f) ' hours'])

EmptyMatrix = zeros(size(r));
C_B = EmptyMatrix+C_B0; % mmol/cm3; fixed BCN concentration
C_P0 = EmptyMatrix; % mmol/cm3; free PEG-4Az concentration
C_PF = EmptyMatrix; % mmol/cm3; newly fixed PEG-4Az concentration
C_X = EmptyMatrix; % mmol/cm3; new XL concentration

C_P = C_P0+C_PF; % mmol/cm3; total new PEG concentration

C_P0_plot = 1000*C_P0; % mM; free PEG-4Az concentration
C_PF_plot = 1000*C_PF; % mM; fixed PEG-4Az concentration
C_P_plot = 1000*C_P;
C_X_plot = 1000*C_X; % mM; new XL concentration (this stiffening step)
C_XT_plot = 1000*(C_XT+C_X); % mM; total new XL concentration over course of entire stiffening

% Heatmap Plot
fig = figure();
fig.Position(3:4) = [1300 650];
hold on
    % free PEG concentration (C_P0)
    ax1 = subplot(3,3,3);
    Plot1 = surf(r,z,C_P0_plot);
    set(Plot1,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(autumn);
    colormap(ax1,cmap)
    title(['Solution Concentration: ' num2str(C_i*1000) ' mM PEG'],'FontSize',14)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'free PEG concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 1.755])
    txt = ['Free PEG'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
    
    % fixed PEG Concentration (C_PF)
    ax2 = subplot(3,3,6);
    Plot2 = surf(r,z,C_PF_plot);
    set(Plot2,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(spring);
    colormap(ax2,cmap)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'fixed PEG concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 1.155])
    txt = ['Fixed PEG'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
    
    % total PEG (C_P)
    ax3 = subplot(3,3,9);
    Plot3 = surf(r,z,C_P_plot);
    set(Plot3,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(copper);
    colormap(ax3,cmap)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'total PEG concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 2.91])
    txt = ['Total PEG'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
    
    % stiffened linkage concentration (C_X)
    ax4 = subplot(3,3,2);
    Plot4 = surf(r,z,C_X_plot);
    set(Plot4,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(gray);
    colormap(ax4,cmap)
    title('Initial Conditions','FontSize',14)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'new linkage concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([-0.05 2.31])
    txt = ['New Linkages (Within Incubation)'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
hold off
frame = getframe(gcf);
writeVideo(writerPEG,frame)

% Heatmap Plot
fig = figure();
fig.Position(3:4) = [800 500];
hold on
    % total linkage concentration (C_XT)
    Plot = surf(r,z,C_XT_plot);
    set(Plot,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    title('Initial Conditions','FontSize',16)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'total stiffened linkage concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 C_XT_max])
    txt = ['Total New Covalent Linkages Over the Course of Stiffening'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
hold off
frame = getframe(gcf);
writeVideo(writerXL,frame)

%% Stiffening Step 2B (PEG-4Azide)
C_inf = C_i; % mmol/cm3 [=] mM/1000;  initial PEG-4Az concentration in solution (250 mL)
disp(['Solution Concentration: ' num2str(C_inf*1000) ' mM PEG'])

PlotTime = [1 : 1 : 12]; % in hours
% Cycle through PlotTime
for b = 1:length(PlotTime)
    time = PlotTime(b);
    t0 = t_f + delta_t; % initial timepoint (in seconds)
    t_f = time*60*60; % plot timepoint in seconds
    t_list = [t0 : delta_t : t_f]; % t steps
    disp(['Timepoint: ' num2str(time) ' hours'])
    
    % Calculate concentrations at each time step
    for a = 1:length(t_list)
        C_B_NEW = C_B-k*delta_t*C_B.*C_P0*u; % new fixed BCN concentration
        C_X_NEW = C_X+k*delta_t*C_B.*C_P0*u; % new XL concentration
        C_PF_NEW = C_PF+k*delta_t*C_B.*C_P0; % new fixed PEG concentration
        
        % Calculate new free PEG concentration
        C_P0_NEW = EmptyMatrix;
        
        % BC at r = R (i = N)
        C_P0_NEW(:,end) = C_inf;
        % BC at z = H (j = M)
        C_P0_NEW(end,:) = C_inf;
        % BC at r = 0 (i = 1), z = 0 (j = 1)
        C_P0_NEW(1,1) = C_P0(1,1) - ((7*D*delta_t*C_P(1,1))/(delta_r^2)+(7*D*delta_t*C_P(1,1))/(2*delta_z^2)+k*delta_t*C_B(1,1)*C_P0(1,1)) + ((8*D*delta_t)/(delta_r^2))*C_P(1,2) - ((D*delta_t)/(delta_r^2))*C_P(1,3) + ((4*D*delta_t)/(delta_z^2))*C_P(2,1) - ((D*delta_t)/(2*delta_z^2))*C_P(3,1);
        % 'BC' (axis of symmetry) at r = 0 (i = 1)
        for j = 2:M-1
            C_P0_NEW(j,1) = C_P0(j,1) - ((7*D*delta_t*C_P(j,1))/(delta_r^2)+(2*D*delta_t*C_P(j,1))/(delta_z^2)+k*delta_t*C_B(j,1)*C_P0(j,1)) + ((8*D*delta_t)/(delta_r^2))*C_P(j,2) - ((D*delta_t)/(delta_r^2))*C_P(j,3) + ((D*delta_t)/(delta_z^2))*C_P(j+1,1) + ((D*delta_t)/(delta_z^2))*C_P(j-1,1);
        end
        % BC at z = 0 (j = 1)
        for i = 2:N-1
            r_i = (i-1)*delta_r;
            C_P0_NEW(1,i) = C_P0(1,i) - ((2*D*delta_t*C_P(1,i))/(delta_r^2)+(7*D*delta_t*C_P(1,i))/(2*delta_z^2)+k*delta_t*C_B(1,i)*C_P0(1,i)) + ((4*D*delta_t)/(delta_z^2))*C_P(2,i) - ((D*delta_t)/(2*delta_z^2))*C_P(3,i) + ((D*delta_t)/(delta_r^2) + (D*delta_t)/(2*r_i*delta_r))*C_P(1,i+1) + ((D*delta_t)/(delta_r^2) - (D*delta_t)/(2*r_i*delta_r))*C_P(1,i-1);
        end
        % internal nodes
        for j = 2:M-1
            for i = 2:N-1
                C_P0_NEW(j,i) = C_P0(j,i) - ((2*D*delta_t*C_P(j,i))/(delta_r^2)+(2*D*delta_t*C_P(j,i))/(delta_z^2)+k*delta_t*C_B(j,i)*C_P0(j,i)) + ((D*delta_t)/(delta_r^2) + (D*delta_t)/(2*r_i*delta_r))*C_P(j,i+1) + ((D*delta_t)/(delta_r^2) - (D*delta_t)/(2*r_i*delta_r))*C_P(j,i-1) + ((D*delta_t)/(delta_z^2))*C_P(j+1,i) + ((D*delta_t)/(delta_z^2))*C_P(j-1,i);
            end
        end
        
        % Update new concentrations
        C_B = C_B_NEW; % mmol/cm3; fixed BCN concentration
        C_P0 = C_P0_NEW; % mmol/cm3; free PEG-4Az concentration
        C_PF = C_PF_NEW; % mmol/cm3; newly fixed PEG-4Az concentration
        C_X = C_X_NEW; % mmol/cm3; new XL concentration
        C_P = C_P0+C_PF; % mmol/cm3; total new PEG concentration
    end
        
    % Update bulk solution concentration
    C_gel_slice = zeros(M-1,N-1);
    for i = 1:N-1
        for j = 1:M-1
            C_gel_slice(j,i) = delta_r*delta_z*((C_P(j,i)+C_P(j,i+1)+C_P(j+1,i)+C_P(j+1,i+1))/4); % mmol/cm3
            % this is the PEG (fixed + free) concentration along a
            % slice of the gel: at each box of delta r and delta z
            % concentration of free PEG + reacted PEG
        end
    end
    C_gel_slice_sum = sum(sum(C_gel_slice)); % equivalent to C_P*R*H (mmol/cm)
    C_gel_full = C_gel_slice_sum*pi*R; % remainder of volume (mmol)
    C_inf = (C_i*V - C_gel_full)/V; % mmol/cm3; new C_inf value for next hour
    
    disp(['New Solution Concentration: ' num2str(C_inf*1000) ' mM PEG'])
    
    % Prepare to plot
    C_P0_plot = 1000*C_P0; % mM; free PEG-4BCN concentration
    C_PF_plot = 1000*C_PF; % mM; fixed PEG-4BCN concentration
    C_P_plot = 1000*C_P; % mM; total PEG-4BCN concentration
    C_X_plot = 1000*C_X; % mM; new XL concentration (this stiffening step)
    C_XT_plot = 1000*(C_XT+C_X); % mM; total new XL concentration over course of entire stiffening
    
    % Heatmap Plot
    fig = figure();
    fig.Position(3:4) = [1300 650];
    hold on
        % free PEG concentration (C_P0)
        ax1 = subplot(3,3,3);
        Plot1 = surf(r,z,C_P0_plot);
        set(Plot1,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(autumn);
        colormap(ax1,cmap)
        title(['Solution Concentration: ' num2str(C_inf*1000) ' mM PEG'],'FontSize',14)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'free PEG concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 1.755])
        txt = ['Free PEG'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off

        % fixed PEG Concentration (C_PF)
        ax2 = subplot(3,3,6);
        Plot2 = surf(r,z,C_PF_plot);
        set(Plot2,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(spring);
        colormap(ax2,cmap)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'fixed PEG concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 1.155])
        txt = ['Fixed PEG'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off

        % total PEG (C_P)
        ax3 = subplot(3,3,9);
        Plot3 = surf(r,z,C_P_plot);
        set(Plot3,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(copper);
        colormap(ax3,cmap)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'total PEG concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 2.91])
        txt = ['Total PEG'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off

        % stiffened linkage concentration (C_X)
        ax4 = subplot(3,3,2);
        Plot4 = surf(r,z,C_X_plot);
        set(Plot4,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(gray);
        colormap(ax4,cmap)
        title([num2str(time) '-Hour Timepoint'],'FontSize',14)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'new linkage concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([-0.05 2.31])
        txt = ['New Linkages (Within Incubation)'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off
    hold off
    frame = getframe(gcf);
    writeVideo(writerPEG,frame)

    % Heatmap Plot
    fig = figure();
    fig.Position(3:4) = [800 500];
    hold on
        % total linkage concentration (C_XT)
        Plot = surf(r,z,C_XT_plot);
        set(Plot,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        title([num2str(time) '-Hour Timepoint'],'FontSize',14)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'total stiffened linkage concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 C_XT_max])
        txt = ['Total New Covalent Linkages Over the Course of Stiffening'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off
    hold off
    frame = getframe(gcf);
    writeVideo(writerXL,frame)
end
C_A = C_X+C_A0; % mmol/cm3; fixed azide concentration (new + leftover from last cycle)
C_XT = C_X+C_XT; % update total linkages over course of stiffening

P0_min = min(min(C_P0))*1000
P0_max = max(max(C_P0))*1000
PF_min = min(min(C_PF))*1000
PF_max = max(max(C_PF))*1000
P_min = min(min(C_P))*1000
P_max = max(max(C_P))*1000
X_min = min(min(C_X))*1000
X_max = max(max(C_X))*1000
XT_min = min(min(C_XT))*1000
XT_max = max(max(C_XT))*1000
%pause(120)

%% Stiffening Step 3A (PEG-4BCN): Initial Conditions

disp(' ');
disp('Stiffening Step 3A: PEG-4BCN');

C_B0 = C_B; % mmol/cm3 [=] mM/1000; fixed BCN functional groups throughout gel
C_A0 = C_A; % mmol/cm3 [=] mM/1000; fixed azide functional groups throughout gel
% Disregard any free PEG remaining in the gel

C_i_BCN = 7.02/1000; % mmol/cm3 [=] mM/1000; initial BCN concentration in solution (250 mL)
C_i = C_i_BCN/4; % mmol/cm3 [=] mM/1000; initial PEG-4BCN concentration in solution (250 mL)

A_max = max(max(C_A0)); % for plotting

% NOTE: treat PEG as though only 2 BCNs will react and the other 2 will be
% fixed to gel
u = 2; % 2 'available' BCN for reaction per PEG
t_f = 0;
disp(['Timepoint: ' num2str(t_f) ' hours'])

EmptyMatrix = zeros(size(r));
C_A = EmptyMatrix+C_A0; % mmol/cm3; fixed azide concentration
C_P0 = EmptyMatrix; % mmol/cm3; free PEG-4BCN concentration
C_PF = EmptyMatrix; % mmol/cm3; newly fixed PEG-4BCN concentration
C_X = EmptyMatrix; % mmol/cm3; new XL concentration

C_P = C_P0+C_PF; % mmol/cm3; total new PEG concentration

C_P0_plot = 1000*C_P0; % mM; free PEG-4Az concentration
C_PF_plot = 1000*C_PF; % mM; fixed PEG-4Az concentration
C_P_plot = 1000*C_P;
C_X_plot = 1000*C_X; % mM; new XL concentration (this stiffening step)
C_XT_plot = 1000*(C_XT+C_X); % mM; total new XL concentration over course of entire stiffening

% Heatmap Plot
fig = figure();
fig.Position(3:4) = [1300 650];
hold on
    % free PEG concentration (C_P0)
    ax1 = subplot(3,3,3);
    Plot1 = surf(r,z,C_P0_plot);
    set(Plot1,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(autumn);
    colormap(ax1,cmap)
    title(['Solution Concentration: ' num2str(C_i*1000) ' mM PEG'],'FontSize',14)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'free PEG concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 1.755])
    txt = ['Free PEG'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
    
    % fixed PEG Concentration (C_PF)
    ax2 = subplot(3,3,6);
    Plot2 = surf(r,z,C_PF_plot);
    set(Plot2,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(spring);
    colormap(ax2,cmap)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'fixed PEG concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 1.155])
    txt = ['Fixed PEG'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
    
    % total PEG (C_P)
    ax3 = subplot(3,3,9);
    Plot3 = surf(r,z,C_P_plot);
    set(Plot3,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(copper);
    colormap(ax3,cmap)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'total PEG concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 2.91])
    txt = ['Total PEG'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
    
    % stiffened linkage concentration (C_X)
    ax4 = subplot(3,3,2);
    Plot4 = surf(r,z,C_X_plot);
    set(Plot4,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(gray);
    colormap(ax4,cmap)
    title('Initial Conditions','FontSize',14)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'new linkage concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([-0.05 2.31])
    txt = ['New Linkages (Within Incubation)'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
hold off
frame = getframe(gcf);
writeVideo(writerPEG,frame)

% Heatmap Plot
fig = figure();
fig.Position(3:4) = [800 500];
hold on
    % total linkage concentration (C_XT)
    Plot = surf(r,z,C_XT_plot);
    set(Plot,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    title('Initial Conditions','FontSize',16)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'total stiffened linkage concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 C_XT_max])
    txt = ['Total New Covalent Linkages Over the Course of Stiffening'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
hold off
frame = getframe(gcf);
writeVideo(writerXL,frame)

%% Stiffening Step 3A (PEG-4BCN)
C_inf = C_i; % mmol/cm3 [=] mM/1000;  initial PEG-4BCN concentration in solution (250 mL)
disp(['Solution Concentration: ' num2str(C_inf*1000) ' mM PEG'])

PlotTime = [1 : 1 : 12]; % in hours
% Cycle through PlotTime
for b = 1:length(PlotTime)
    time = PlotTime(b);
    t0 = t_f + delta_t; % initial timepoint (in seconds)
    t_f = time*60*60; % plot timepoint in seconds
    t_list = [t0 : delta_t : t_f]; % t steps
    disp(['Timepoint: ' num2str(time) ' hours'])
    
    % Calculate concentrations at each time step
    for a = 1:length(t_list)
        C_A_NEW = C_A-k*delta_t*C_A.*C_P0*u; % new fixed azide concentration
        C_X_NEW = C_X+k*delta_t*C_A.*C_P0*u; % new XL concentration
        C_PF_NEW = C_PF+k*delta_t*C_A.*C_P0; % new fixed PEG concentration
        
        % Calculate new free PEG concentration
        C_P0_NEW = EmptyMatrix;
        
        % BC at r = R (i = N)
        C_P0_NEW(:,end) = C_inf;
        % BC at z = H (j = M)
        C_P0_NEW(end,:) = C_inf;
        % BC at r = 0 (i = 1), z = 0 (j = 1)
        C_P0_NEW(1,1) = C_P0(1,1) - ((7*D*delta_t*C_P(1,1))/(delta_r^2)+(7*D*delta_t*C_P(1,1))/(2*delta_z^2)+k*delta_t*C_A(1,1)*C_P0(1,1)) + ((8*D*delta_t)/(delta_r^2))*C_P(1,2) - ((D*delta_t)/(delta_r^2))*C_P(1,3) + ((4*D*delta_t)/(delta_z^2))*C_P(2,1) - ((D*delta_t)/(2*delta_z^2))*C_P(3,1);
        % 'BC' (axis of symmetry) at r = 0 (i = 1)
        for j = 2:M-1
            C_P0_NEW(j,1) = C_P0(j,1) - ((7*D*delta_t*C_P(j,1))/(delta_r^2)+(2*D*delta_t*C_P(j,1))/(delta_z^2)+k*delta_t*C_A(j,1)*C_P0(j,1)) + ((8*D*delta_t)/(delta_r^2))*C_P(j,2) - ((D*delta_t)/(delta_r^2))*C_P(j,3) + ((D*delta_t)/(delta_z^2))*C_P(j+1,1) + ((D*delta_t)/(delta_z^2))*C_P(j-1,1);
        end
        % BC at z = 0 (j = 1)
        for i = 2:N-1
            r_i = (i-1)*delta_r;
            C_P0_NEW(1,i) = C_P0(1,i) - ((2*D*delta_t*C_P(1,i))/(delta_r^2)+(7*D*delta_t*C_P(1,i))/(2*delta_z^2)+k*delta_t*C_A(1,i)*C_P0(1,i)) + ((4*D*delta_t)/(delta_z^2))*C_P(2,i) - ((D*delta_t)/(2*delta_z^2))*C_P(3,i) + ((D*delta_t)/(delta_r^2) + (D*delta_t)/(2*r_i*delta_r))*C_P(1,i+1) + ((D*delta_t)/(delta_r^2) - (D*delta_t)/(2*r_i*delta_r))*C_P(1,i-1);
        end
        % internal nodes
        for j = 2:M-1
            for i = 2:N-1
                C_P0_NEW(j,i) = C_P0(j,i) - ((2*D*delta_t*C_P(j,i))/(delta_r^2)+(2*D*delta_t*C_P(j,i))/(delta_z^2)+k*delta_t*C_A(j,i)*C_P0(j,i)) + ((D*delta_t)/(delta_r^2) + (D*delta_t)/(2*r_i*delta_r))*C_P(j,i+1) + ((D*delta_t)/(delta_r^2) - (D*delta_t)/(2*r_i*delta_r))*C_P(j,i-1) + ((D*delta_t)/(delta_z^2))*C_P(j+1,i) + ((D*delta_t)/(delta_z^2))*C_P(j-1,i);
            end
        end
        
        % Update new concentrations
        C_A = C_A_NEW; % mmol/cm3; fixed azide concentration
        C_P0 = C_P0_NEW; % mmol/cm3; free PEG-4BCN concentration
        C_PF = C_PF_NEW; % mmol/cm3; newly fixed PEG-4BCN concentration
        C_X = C_X_NEW; % mmol/cm3; new XL concentration
        C_P = C_P0+C_PF; % mmol/cm3; total new PEG concentration
    end
        
    % Update bulk solution concentration if it's been an hour
    C_gel_slice = zeros(M-1,N-1);
    for i = 1:N-1
        for j = 1:M-1
            C_gel_slice(j,i) = delta_r*delta_z*((C_P(j,i)+C_P(j,i+1)+C_P(j+1,i)+C_P(j+1,i+1))/4); % mmol/cm3
            % this is the PEG (fixed + free) concentration along a
            % slice of the gel: at each box of delta r and delta z
            % concentration of free PEG + reacted PEG
        end
    end
    C_gel_slice_sum = sum(sum(C_gel_slice)); % equivalent to C_P*R*H (mmol/cm)
    C_gel_full = C_gel_slice_sum*pi*R; % remainder of volume (mmol)
    C_inf = (C_i*V - C_gel_full)/V; % mmol/cm3; new C_inf value for next hour

    disp(['New Solution Concentration: ' num2str(C_inf*1000) ' mM PEG'])
    
    % Prepare to plot
    C_P0_plot = 1000*C_P0; % mM; free PEG-4BCN concentration
    C_PF_plot = 1000*C_PF; % mM; fixed PEG-4BCN concentration
    C_P_plot = 1000*C_P; % mM; total PEG-4BCN concentration
    C_X_plot = 1000*C_X; % mM; new XL concentration (this stiffening step)
    C_XT_plot = 1000*(C_XT+C_X); % mM; total new XL concentration over course of entire stiffening
    
    % Heatmap Plot
    fig = figure();
    fig.Position(3:4) = [1300 650];
    hold on
        % free PEG concentration (C_P0)
        ax1 = subplot(3,3,3);
        Plot1 = surf(r,z,C_P0_plot);
        set(Plot1,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(autumn);
        colormap(ax1,cmap)
        title(['Solution Concentration: ' num2str(C_inf*1000) ' mM PEG'],'FontSize',14)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'free PEG concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 1.755])
        txt = ['Free PEG'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off

        % fixed PEG Concentration (C_PF)
        ax2 = subplot(3,3,6);
        Plot2 = surf(r,z,C_PF_plot);
        set(Plot2,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(spring);
        colormap(ax2,cmap)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'fixed PEG concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 1.155])
        txt = ['Fixed PEG'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off

        % total PEG (C_P)
        ax3 = subplot(3,3,9);
        Plot3 = surf(r,z,C_P_plot);
        set(Plot3,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(copper);
        colormap(ax3,cmap)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'total PEG concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 2.91])
        txt = ['Total PEG'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off

        % stiffened linkage concentration (C_X)
        ax4 = subplot(3,3,2);
        Plot4 = surf(r,z,C_X_plot);
        set(Plot4,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(gray);
        colormap(ax4,cmap)
        title([num2str(time) '-Hour Timepoint'],'FontSize',14)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'new linkage concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([-0.05 2.31])
        txt = ['New Linkages (Within Incubation)'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off
    hold off
    frame = getframe(gcf);
    writeVideo(writerPEG,frame)

    % Heatmap Plot
    fig = figure();
    fig.Position(3:4) = [800 500];
    hold on
        % total linkage concentration (C_XT)
        Plot = surf(r,z,C_XT_plot);
        set(Plot,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        title([num2str(time) '-Hour Timepoint'],'FontSize',14)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'total stiffened linkage concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 C_XT_max])
        txt = ['Total New Covalent Linkages Over the Course of Stiffening'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off
    hold off
    frame = getframe(gcf);
    writeVideo(writerXL,frame)
end
C_B = C_X+C_B0; % mmol/cm3; fixed BCN concentration
C_XT = C_X+C_XT; % update total linkages over course of stiffening

P0_min = min(min(C_P0))*1000
P0_max = max(max(C_P0))*1000
PF_min = min(min(C_PF))*1000
PF_max = max(max(C_PF))*1000
P_min = min(min(C_P))*1000
P_max = max(max(C_P))*1000
X_min = min(min(C_X))*1000
X_max = max(max(C_X))*1000
XT_min = min(min(C_XT))*1000
XT_max = max(max(C_XT))*1000
%pause(120)

%% Stiffening Step 3B (PEG-4Azide): Initial Conditions

disp(' ');
disp('Stiffening Step 3B: PEG-4Azide');

C_B0 = C_B; % mmol/cm3 [=] mM/1000; fixed BCN functional groups throughout gel
C_A0 = C_A; % mmol/cm3 [=] mM/1000; fixed azide functional groups throughout gel
% Disregard any free PEG remaining in the gel

C_i_Az = 3.51/1000; % mmol/cm3 [=] mM/1000; initial azide concentration in solution (250 mL)
C_i = C_i_Az/4; % mmol/cm3 [=] mM/1000; initial PEG-4Az concentration in solution (250 mL)

B_max = max(max(C_B0)); % for plotting

% NOTE: treat PEG as though only 2 azides will react and the other 2 will
% be fixed to gel
u = 4; % 4 'available' azide for reaction per PEG
t_f = 0;
disp(['Timepoint: ' num2str(t_f) ' hours'])

EmptyMatrix = zeros(size(r));
C_B = EmptyMatrix+C_B0; % mmol/cm3; fixed BCN concentration
C_P0 = EmptyMatrix; % mmol/cm3; free PEG-4Az concentration
C_PF = EmptyMatrix; % mmol/cm3; newly fixed PEG-4Az concentration
C_X = EmptyMatrix; % mmol/cm3; new XL concentration

C_P = C_P0+C_PF; % mmol/cm3; total new PEG concentration

C_P0_plot = 1000*C_P0; % mM; free PEG-4Az concentration
C_PF_plot = 1000*C_PF; % mM; fixed PEG-4Az concentration
C_P_plot = 1000*C_P;
C_X_plot = 1000*C_X; % mM; new XL concentration (this stiffening step)
C_XT_plot = 1000*(C_XT+C_X); % mM; total new XL concentration over course of entire stiffening

% Heatmap Plot
fig = figure();
fig.Position(3:4) = [1300 650];
hold on
    % free PEG concentration (C_P0)
    ax1 = subplot(3,3,3);
    Plot1 = surf(r,z,C_P0_plot);
    set(Plot1,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(autumn);
    colormap(ax1,cmap)
    title(['Solution Concentration: ' num2str(C_i*1000) ' mM PEG'],'FontSize',14)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'free PEG concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 1.755])
    txt = ['Free PEG'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
    
    % fixed PEG Concentration (C_PF)
    ax2 = subplot(3,3,6);
    Plot2 = surf(r,z,C_PF_plot);
    set(Plot2,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(spring);
    colormap(ax2,cmap)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'fixed PEG concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 1.155])
    txt = ['Fixed PEG'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
    
    % total PEG (C_P)
    ax3 = subplot(3,3,9);
    Plot3 = surf(r,z,C_P_plot);
    set(Plot3,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(copper);
    colormap(ax3,cmap)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'total PEG concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 2.91])
    txt = ['Total PEG'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
    
    % stiffened linkage concentration (C_X)
    ax4 = subplot(3,3,2);
    Plot4 = surf(r,z,C_X_plot);
    set(Plot4,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    cmap = flip(gray);
    colormap(ax4,cmap)
    title('Initial Conditions','FontSize',14)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'new linkage concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([-0.05 2.31])
    txt = ['New Linkages (Within Incubation)'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
hold off
frame = getframe(gcf);
writeVideo(writerPEG,frame)

% Heatmap Plot
fig = figure();
fig.Position(3:4) = [800 500];
hold on
    % total linkage concentration (C_XT)
    Plot = surf(r,z,C_XT_plot);
    set(Plot,'edgecolor','none')
    daspect([1 1 1])
    view(2)
    title('Initial Conditions','FontSize',16)
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'total stiffened linkage concentration [mM]','Rotation',-270,'FontSize',12);
    xlim([0 0.47])
    ylim([0 0.25])
    caxis([0 C_XT_max])
    txt = ['Total New Covalent Linkages Over the Course of Stiffening'];
    text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
    grid off
hold off
frame = getframe(gcf);
writeVideo(writerXL,frame)

%% Stiffening Step 3B (PEG-4Azide)
C_inf = C_i; % mmol/cm3 [=] mM/1000;  initial PEG-4Az concentration in solution (250 mL)
disp(['Solution Concentration: ' num2str(C_inf*1000) ' mM PEG'])

PlotTime = [1 : 1 : 12]; % in hours
% Cycle through PlotTime
for b = 1:length(PlotTime)
    time = PlotTime(b);
    t0 = t_f + delta_t; % initial timepoint (in seconds)
    t_f = time*60*60; % plot timepoint in seconds
    t_list = [t0 : delta_t : t_f]; % t steps
    disp(['Timepoint: ' num2str(time) ' hours'])
    
    % Calculate concentrations at each time step
    for a = 1:length(t_list)
        C_B_NEW = C_B-k*delta_t*C_B.*C_P0*u; % new fixed BCN concentration
        C_X_NEW = C_X+k*delta_t*C_B.*C_P0*u; % new XL concentration
        C_PF_NEW = C_PF+k*delta_t*C_B.*C_P0; % new fixed PEG concentration
        
        % Calculate new free PEG concentration
        C_P0_NEW = EmptyMatrix;
        
        % BC at r = R (i = N)
        C_P0_NEW(:,end) = C_inf;
        % BC at z = H (j = M)
        C_P0_NEW(end,:) = C_inf;
        % BC at r = 0 (i = 1), z = 0 (j = 1)
        C_P0_NEW(1,1) = C_P0(1,1) - ((7*D*delta_t*C_P(1,1))/(delta_r^2)+(7*D*delta_t*C_P(1,1))/(2*delta_z^2)+k*delta_t*C_B(1,1)*C_P0(1,1)) + ((8*D*delta_t)/(delta_r^2))*C_P(1,2) - ((D*delta_t)/(delta_r^2))*C_P(1,3) + ((4*D*delta_t)/(delta_z^2))*C_P(2,1) - ((D*delta_t)/(2*delta_z^2))*C_P(3,1);
        % 'BC' (axis of symmetry) at r = 0 (i = 1)
        for j = 2:M-1
            C_P0_NEW(j,1) = C_P0(j,1) - ((7*D*delta_t*C_P(j,1))/(delta_r^2)+(2*D*delta_t*C_P(j,1))/(delta_z^2)+k*delta_t*C_B(j,1)*C_P0(j,1)) + ((8*D*delta_t)/(delta_r^2))*C_P(j,2) - ((D*delta_t)/(delta_r^2))*C_P(j,3) + ((D*delta_t)/(delta_z^2))*C_P(j+1,1) + ((D*delta_t)/(delta_z^2))*C_P(j-1,1);
        end
        % BC at z = 0 (j = 1)
        for i = 2:N-1
            r_i = (i-1)*delta_r;
            C_P0_NEW(1,i) = C_P0(1,i) - ((2*D*delta_t*C_P(1,i))/(delta_r^2)+(7*D*delta_t*C_P(1,i))/(2*delta_z^2)+k*delta_t*C_B(1,i)*C_P0(1,i)) + ((4*D*delta_t)/(delta_z^2))*C_P(2,i) - ((D*delta_t)/(2*delta_z^2))*C_P(3,i) + ((D*delta_t)/(delta_r^2) + (D*delta_t)/(2*r_i*delta_r))*C_P(1,i+1) + ((D*delta_t)/(delta_r^2) - (D*delta_t)/(2*r_i*delta_r))*C_P(1,i-1);
        end
        % internal nodes
        for j = 2:M-1
            for i = 2:N-1
                C_P0_NEW(j,i) = C_P0(j,i) - ((2*D*delta_t*C_P(j,i))/(delta_r^2)+(2*D*delta_t*C_P(j,i))/(delta_z^2)+k*delta_t*C_B(j,i)*C_P0(j,i)) + ((D*delta_t)/(delta_r^2) + (D*delta_t)/(2*r_i*delta_r))*C_P(j,i+1) + ((D*delta_t)/(delta_r^2) - (D*delta_t)/(2*r_i*delta_r))*C_P(j,i-1) + ((D*delta_t)/(delta_z^2))*C_P(j+1,i) + ((D*delta_t)/(delta_z^2))*C_P(j-1,i);
            end
        end
        
        % Update new concentrations
        C_B = C_B_NEW; % mmol/cm3; fixed BCN concentration
        C_P0 = C_P0_NEW; % mmol/cm3; free PEG-4Az concentration
        C_PF = C_PF_NEW; % mmol/cm3; newly fixed PEG-4Az concentration
        C_X = C_X_NEW; % mmol/cm3; new XL concentration
        C_P = C_P0+C_PF; % mmol/cm3; total new PEG concentration
    end
        
    % Update bulk solution concentration
    C_gel_slice = zeros(M-1,N-1);
    for i = 1:N-1
        for j = 1:M-1
            C_gel_slice(j,i) = delta_r*delta_z*((C_P(j,i)+C_P(j,i+1)+C_P(j+1,i)+C_P(j+1,i+1))/4); % mmol/cm3
            % this is the PEG (fixed + free) concentration along a
            % slice of the gel: at each box of delta r and delta z
            % concentration of free PEG + reacted PEG
        end
    end
    C_gel_slice_sum = sum(sum(C_gel_slice)); % equivalent to C_P*R*H (mmol/cm)
    C_gel_full = C_gel_slice_sum*pi*R; % remainder of volume (mmol)
    C_inf = (C_i*V - C_gel_full)/V; % mmol/cm3; new C_inf value for next hour
        
    disp(['New Solution Concentration: ' num2str(C_inf*1000) ' mM PEG'])
    
    % Prepare to plot
    C_P0_plot = 1000*C_P0; % mM; free PEG-4BCN concentration
    C_PF_plot = 1000*C_PF; % mM; fixed PEG-4BCN concentration
    C_P_plot = 1000*C_P; % mM; total PEG-4BCN concentration
    C_X_plot = 1000*C_X; % mM; new XL concentration (this stiffening step)
    C_XT_plot = 1000*(C_XT+C_X); % mM; total new XL concentration over course of entire stiffening
    
    % Heatmap Plot
    fig = figure();
    fig.Position(3:4) = [1300 650];
    hold on
        % free PEG concentration (C_P0)
        ax1 = subplot(3,3,3);
        Plot1 = surf(r,z,C_P0_plot);
        set(Plot1,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(autumn);
        colormap(ax1,cmap)
        title(['Solution Concentration: ' num2str(C_inf*1000) ' mM PEG'],'FontSize',14)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'free PEG concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 1.755])
        txt = ['Free PEG'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off

        % fixed PEG Concentration (C_PF)
        ax2 = subplot(3,3,6);
        Plot2 = surf(r,z,C_PF_plot);
        set(Plot2,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(spring);
        colormap(ax2,cmap)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'fixed PEG concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 1.155])
        txt = ['Fixed PEG'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off

        % total PEG (C_P)
        ax3 = subplot(3,3,9);
        Plot3 = surf(r,z,C_P_plot);
        set(Plot3,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(copper);
        colormap(ax3,cmap)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'total PEG concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 2.91])
        txt = ['Total PEG'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off

        % stiffened linkage concentration (C_X)
        ax4 = subplot(3,3,2);
        Plot4 = surf(r,z,C_X_plot);
        set(Plot4,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        cmap = flip(gray);
        colormap(ax4,cmap)
        title([num2str(time) '-Hour Timepoint'],'FontSize',14)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'new linkage concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([-0.05 2.31])
        txt = ['New Linkages (Within Incubation)'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off
    hold off
    frame = getframe(gcf);
    writeVideo(writerPEG,frame)

    % Heatmap Plot
    fig = figure();
    fig.Position(3:4) = [800 500];
    hold on
        % total linkage concentration (C_XT)
        Plot = surf(r,z,C_XT_plot);
        set(Plot,'edgecolor','none')
        daspect([1 1 1])
        view(2)
        title([num2str(time) '-Hour Timepoint'],'FontSize',14)
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'total stiffened linkage concentration [mM]','Rotation',-270,'FontSize',12);
        xlim([0 0.47])
        ylim([0 0.25])
        caxis([0 C_XT_max])
        txt = ['Total New Covalent Linkages Over the Course of Stiffening'];
        text(0.02,0.225,txt,'FontSize',12,'FontWeight','bold')
        grid off
    hold off
    frame = getframe(gcf);
    writeVideo(writerXL,frame)
end
close(writerPEG)
close(writerXL)
C_A = C_X-4*C_PF+C_A0; % mmol/cm3; fixed azide concentration (new + leftover from last cycle)
C_XT = C_X+C_XT; % update total linkages over course of stiffening

P0_min = min(min(C_P0))*1000
P0_max = max(max(C_P0))*1000
PF_min = min(min(C_PF))*1000
PF_max = max(max(C_PF))*1000
P_min = min(min(C_P))*1000
P_max = max(max(C_P))*1000
X_min = min(min(C_X))*1000
X_max = max(max(C_X))*1000
XT_min = min(min(C_XT))*1000
XT_max = max(max(C_XT))*1000